## Libraries
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
# library(purrr)
# library(httr)
# library(magrittr)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
library(rgdal)
## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files:
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 7.0.0, March 1st, 2020, [PJ_VERSION: 700]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.4-1
library(lidR)
library(functional)
library(RColorBrewer)
library(glue)
##
## Attaching package: 'glue'
## The following object is masked from 'package:raster':
##
## trim
## The following object is masked from 'package:dplyr':
##
## collapse
library(furrr)
## Loading required package: future
##
## Attaching package: 'future'
## The following object is masked from 'package:raster':
##
## values
# future::plan(multiprocess)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
crs <- st_crs(32650)
crs$proj4string
## [1] "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
years <- c(2013,2014)
ws <-c(15)
dir_data <-"../Data/lidar"
dir_chm <- paste(dir_data,"raster",sep="/")
dir_treetops <- paste(dir_data,"vector","treetops",years,sep="/")
dir_crowns <- paste(dir_data,"raster","crowns",years,sep = "/")
walk(dir_crowns,dir.create)
## Warning in .f(.x[[i]], ...): '../Data/lidar/raster/crowns/2013' already exists
## Warning in .f(.x[[i]], ...): '../Data/lidar/raster/crowns/2014' already exists
filename_crowns_index <- paste(dir_data,"raster","crowns","index.txt",sep="/")
filename_crowns_errors <- paste(dir_data,"raster","crowns","errors.txt",sep="/")
## Read rasters
chm_filenames <- paste0(dir_chm,"/raster",years,".tif")
chms <- map(chm_filenames,~raster(.,crs=crs$proj4string))
read_treetops <- function(year,ws){
filename <-paste0("treetops_lmf_ws",ws,".json")
filename_full <- paste(dir_data,"vector","treetops",year,filename,sep = "/")
readOGR(filename_full,p4s = crs$proj4string)
}
treetops <-
map(years,~read_treetops(.,ws))
## OGR data source with driver: GeoJSON
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/treetops/2013/treetops_lmf_ws15.json", layer: "2013_ws15"
## with 3172 features
## It has 2 fields
## Warning in readOGR(filename_full, p4s = crs$proj4string): p4s= argument given as: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
## and read as: +proj=longlat +datum=WGS84 +no_defs
## read string overridden by given p4s= argument value
## OGR data source with driver: GeoJSON
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/treetops/2014/treetops_lmf_ws15.json", layer: "2014_ws15"
## with 3241 features
## It has 2 fields
## Warning in readOGR(filename_full, p4s = crs$proj4string): p4s= argument given as: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
## and read as: +proj=longlat +datum=WGS84 +no_defs
## read string overridden by given p4s= argument value
# # subset tall trees
treetops <-
treetops %>%
# map(~subset(.,.$Z > 60))%>%
map(st_as_sf)%>%
{.}
treetops
## [[1]]
## Simple feature collection with 3172 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 587149.8 ymin: 547214.2 xmax: 588598.8 ymax: 548713.2
## CRS: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
## treeID Z geometry
## 1 1 50.83929 POINT (587182.8 548713.2)
## 2 2 66.56500 POINT (587307.8 548713.2)
## 3 3 65.89000 POINT (587513.2 548713.2)
## 4 4 46.18609 POINT (587549.2 548713.2)
## 5 5 66.71000 POINT (587596.2 548713.2)
## 6 6 62.41000 POINT (587615.8 548713.2)
## 7 7 45.32000 POINT (587692.8 548713.2)
## 8 8 60.35218 POINT (587780.8 548713.2)
## 9 9 46.38237 POINT (587913.2 548713.2)
## 10 10 51.00000 POINT (587933.2 548713.2)
##
## [[2]]
## Simple feature collection with 3241 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 587149.2 ymin: 547213.2 xmax: 588599.2 ymax: 548713.8
## CRS: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
## treeID Z geometry
## 1 1 50.31984 POINT (587182.8 548713.8)
## 2 2 51.99703 POINT (587211.2 548713.8)
## 3 3 46.25211 POINT (587250.8 548713.8)
## 4 4 61.81390 POINT (587394.2 548713.8)
## 5 5 54.84300 POINT (587459.2 548713.8)
## 6 6 64.12300 POINT (587478.8 548713.8)
## 7 7 64.88100 POINT (587514.2 548713.8)
## 8 8 63.59278 POINT (587601.8 548713.8)
## 9 9 37.66848 POINT (587678.8 548713.8)
## 10 10 37.49701 POINT (587692.8 548713.8)
# plot(chms[[1]])
# plot(treetops[[1]]$geometry,add=T)
# plot(chms[[2]])
# plot(treetops[[2]]$geometry,add=T)
plot(chms[[1]])
plot(treetops[[1]]$geometry,add=T)
plot(treetops[[2]]$geometry,col="red",pch="+",add=T)
dm <- st_distance(treetops[[1]],treetops[[2]])
dim(dm)
## [1] 3172 3241
a <-
map(array_branch(dm,1), sort) %>%
as_vector() %>%
array(dim=rev(dim(dm))) %>%
t() %>%
# order()%>%
{.}
b <-a[order(a[,1]),]
plot(b[,1],pch=".")
points(b[,2],pch='.',col='red')
title(paste("Distance to two closest trees, ordered by 1st, size",ws))
b <-a[order(a[,2]),]
plot(b[,1],pch=".")
points(b[,2],pch='.',col='red')
title(paste("Distance to two closest trees, ordered by 2nd, size",ws))
